%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; clc;
date = {'23062020'}
%% manages extraction of data on employment bins
dir = "ThirdChapter_revision";
input = ["/Volumes/Transcend/AWFP/",dir,"/data/CalibrationInput/"];
input = join(input,"");
output= ["/Volumes/Transcend/AWFP/",dir,"/data/CalibrationOutput/"];
output= join(output,"");

cd(input)
options = optimset('Display','off');

E       = exp(csvread('employment_avg_estimation.csv'));            % average employment after s.a. and detrending during 1990-2004
N       = exp(csvread('count_avg_estimation.csv'));                 % average number of firms after s.a. and detrending during 1990-2004
d       = csvread('delta_avg_estimation.csv');                      % Pareto tail from 1990 to 2004
share   = csvread('meanshare_avg_estimation.csv');                  % share of the largest firm, average value during 1990-2004
Ed      = exp(csvread('employment_ts_estimation.csv'));
%%
param = []
Alpha = linspace(.8,.85,3);
Psi   = [1.1125,1.1250,1.1375,1.1500,1.1625];
Gamma = [1];

auxg = repmat(Gamma',size(Alpha,2)*size(Psi,2),1);
auxa = repmat(Alpha,size(Gamma,2)*size(Psi,2),1);
auxa = reshape(auxa,size(Gamma,2)*size(Alpha,2)*size(Psi,2),1);
auxp = repmat(Psi,size(Gamma,2),1);
auxp = reshape(auxp,size(Psi,2)*size(Gamma,2),1);
auxp = repmat(auxp,size(Alpha,2),1);
auxPar = [auxa,auxp,auxg];

niter = size(Gamma,2)*size(Alpha,2)*size(Psi,2);

for y = 1:niter
cd(output)

alpha = auxPar(y,1);
psi   = auxPar(y,2);
gamma = auxPar(y,3);

param(y,:) = [alpha,psi,gamma];
delta = d./(1-alpha);

% setting grid to search for value of S:
Sint = linspace(2,50,49)';                  % prevents from selecting S = 1
Fint = getS(Sint,N,share,delta,alpha,psi);

for i = 1:72
   [V(1,i),S(1,i)] = min(Fint(:,i)); 
end
S = S + 1;

ms = mean(S);
disp(sprintf('mean number of steps is %d',ms))

Ad1= getA(Ed,N,alpha,gamma);                         % get productivity given Ld data        
Am = getAm(N,S,alpha,psi,delta);                     % get model productivity given Pareto tail data

Ad   = Ad1.*(Am./mean(Ad1,1));                       % rescale down A-dyn to get proper shares
int = zeros(size(Ed,1),max(S),size(Ed,2));

for t = 1:size(Ed,1)
for i = 1:size(Ed,2)  
   for s = 1:S(1,i)
      int(t,s,i) = round(psi^(s/(1-alpha)).*(Ed(t,i)/Ad(t,i)));
   end
end
end

for i = 1:size(Ed,2)
    aux = int(:,:,i);
    
    stub = sprintf('/data/CalibrationOutput/Intervals_%s/',date{1});
    path = ['/Volumes/Transcend/AWFP/',dir,stub];
    path = join(path,"");
    cd(path)
    filename = sprintf('interval%dtrial%d.csv',i,y);
    csvwrite(filename,aux);
end

SS(y,:) = S;
VV(y,:) = V;
clear gamma alpha psi
end
%% saving
cd(output)
param = [linspace(1,niter,niter)',param];
filename = sprintf('param%s.csv',date{1})
csvwrite(filename,param)
filename = sprintf('Sint%s.csv',date{1})
csvwrite(filename,SS)
filename = sprintf('param%s.mat',date{1})
save(filename,'param')
filename = sprintf('Sint%s.mat',date{1})
save(filename,'SS')
